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Abstract. Some time ago, Cuniberti et al have proposed a novel method for analytically continuing thermal 
imaginary-time correlators to real time, which requires no model input and should be applicable with finite- 
precision data as well. Given that these assertions go against common wisdom, we report on a naive test 
of the method with an idealized example. We do encounter two problems, which we spell out in detail; this 
implies that systematic errors are difficult to quantify. On a more positive note, the method is simple to 
implement and allows for an empirical recipe by which a reasonable qualitative estimate for some transport 
coefficient may be obtained, if statistical errors of an ultraviolet-subtracted imaginary-time measurement 
can be reduced to roughly below the per mille level. 
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It is a perennial problem that a reliable study of strong 
interactions at temperatures around a few hundred MeV 
requires non-perturbative methods, yet standard Monte 
Carlo techniques only work in Euclidean signature. This 
implies that many of the physically most interesting ob- 
servables, such as transport coefficients or particle pro- 
duction rates (for reviews see e.g. refs. dJ[2]), which are 
inherently Minkowskian in nature, are difficult to address. 

A few years ago, a possible solution to this problem 
was proposed [3], through the use of a "Maximum En- 
tropy Method" 0]. Unfortunately, despite many imple- 
mentations and various related attempts (see e.g. refs. [5]- 
|llj ; similar discussions continue also in the context of 
condensed matter physics problems, see e.g. ref. |12j and 
references therein), it remains a problem that there ap- 
pears to be uncontrolled dependence on model input in 
these results. Therefore, it is perhaps worthwhile to look 
for alternatives as well. 

Some time ago, Cuniberti et al [13] put forward a con- 
crete suggestion which could help in this respect. An al- 
gorithm was provided which was shown to yield a correct 
analytic continuation in a specific limit; furthermore, it 
was suggested that the method might work even if there 
is a finite amount of data and the data points have error 
bars. Despite these attractive features, we are not aware of 
a previous numerical implementation of the algorithm. For 
the record, we wish to report one in this short note, even 
if quantitative success is somewhat marginal. Our hope is 
that some of our theoretical remarks are nevertheless of 
interest and that, on the qualitative level, the algorithm 
might turn out to be a useful addition to the existing tool 
kit. 



2 Basic idea 

The general philosophy of the method can be summarized 
as follows. We consider a periodic ("bosonic") Euclidean 
correlator, Q(r, •), with Q(r + k(3, •) = Q(t, •) where k £ TL, 
/3 denotes the inverse temperature, and • denotes sup- 
pressed variables such as spatial momentum. The correla- 
tor should be analytic everywhere except at Re r = 0+fc/3; 
there it should still be continuous. We also assume the 
further property Q((3 — t, •) = G(t, •), < r < /3, which 
was not imposed in ref. [13] but is satisfied by typical 
gauge-invariant current-current correlators measured in 
lattice QCD. Given the finiteness of /3, the Fourier repre- 
sentation involves a discrete set of Matsubara frequencies, 

g(w n , •) = /,f dr e l ^ T G(T, •), CJ n = 2ttuT, nE%,T = ^~ x . 
Making use of a general hypothesis about the asymptotic 
behaviour of various correlators in the Minkowskian time 
domain, namely that they show at most powerlike growth 
at infinity (physically relevant correlators are actually ex- 
pected to vanish at infinity, see e.g. ref. [2]), a unique 
analytic continuation to a part of the complex plane can 
be shown to exist and can be constructed explicitly. 

More concretely^ the basic quantities defined 
in ref. [T3] are the one-sided sum G + (t,-) = 
TJ2u >o Gfan, •)e _tw * ,T , which is analytic for Imr < 
but has cuts in the upper half-plane; its discontinuity 
across the cut starting at the origin, J + (i, •) = i[£/ + (e + 
it, •) — G + {— e + it, •)], t > 0, e = + , which equals the re- 
tarded real-time correlator, lZ(t, •); as well as its Laplace 
transform J + ((, •) = / dt e _< »* J + (t, ■), which is analytic 
for Ke( > 0. For ( = w n , J + (C, •) reduces to the Fourier 



1 The following relations, though simple to state, may not 
be entirely obvious at first sight; mathematical proofs can be 
found in ref. [13] and references therein. 
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Fig. 1. The complex planes, basic functions, and analytic structures as discussed in the text. 



components G(ui n , •), and therefore constitutes the desired 
analytic continuation to a complex half-plane. The value of 
J + (C, •) along the axis ( = e-jw,w6]R, yields the Fourier 
transform of the retarded correlator, 1Z(uj, •), whose imag- 
inary part in turn equals the spectral function. The basic 
analytic structure is illustrated in fig. [1] 



3 Algorithm 

To implement the analytic continuation, the idea of 
ref. [13] is to expand J + (£, •) with the help of Pollaczek 
polynomials (of the type defined on an infinite interval); 
the retarded correlator J + {t, •) is in turn expressed as 
a linear combination of Laguerre polynomials, with ar- 
gument 2e~ 27TtT . More explicitly, taking ref. [13] at face 
value, the steps are as follows: 

(i) Compute the Fourier modes, G(Lu n , •)• Due to Q((3 — 
t, •) = Q(t, •) S R, the Fourier modes are real and 
even in uj n 

(ii) Construct the coefficients 



^— ' n! 

n=0 



12^1 



K, 71 ■ 



i;i;2) , 

0,1,2,.... (f) 



We have re-expressed the Pollaczek polynomials of 
ref. [T3] through the hypergeometric function 2-F1GI 
Note that the Matsubara zero-mode does not con- 
tribute in eq. ([1]). 

(iii) According to ref. [13 , ae decreases with £, such that 

S?!lo \ ae 1 2 i s finite. 

(iv) Defining t = 2irtT, the retarded real-time correlator 
can now be obtained as 



J+(V) 



ai Li(2e 



(2) 



2 More precisely: Pi(f]) 
ir/;l;2), cf. ref. [IS]. 



where are the Laguerre polynomials. Given that 
Le(0) — 1, the asymptotic value is J + (oo,-) = 
Sfco a ti which should vanish in physically meaningful 
examples [Tl] . 

(v) Given J + , the spectral correlator reads pit) = 
± [6(t)J + (t)-6(-t)J+(-t)] , where J+ denotes a com- 
plex conjugate. Noting that J + as produced by eq. ([2]) 
is real under our assumptions, the spectral function 
can finally be obtained as 



dt sin(tl)<) J + (t, •) , Q 



2ttT 



(3) 



4 Practical implementation 

To apply the previous algorithm to a situation where only 
a finite number of data points are available, we adopt the 
following steps, with the same numbering as in section 3. 

(i) We assume the interval < t < (3 to be evenly divided 
into N parts. The Fourier modes are constructed by 
a discrete transformation; we leave out the "contact 
point" (r = or r = (3), given its possible inaccessi- 
bility to practical measurement. Then, 



S(w n , •) 



N-l 

fe=i 



n = 0, . . . , N - 1 . 

We stress again that, thanks to 9(j3—T, •) = Q{t, 
for < r < /3, Q(u n ,-) £ EL 



(4) 
G R 



(ii) Defining the coefficients </>^„ = 2F\{—i,n + 

1; 1; 2), n > 0, which can be constructed from the re- 
currence relation 



-i = 0; 4>e t a = 1 ( 

{21 + l)^ e ,n-l + 0£, n _ 2 



rr 



, n > 1 , (5) 



the can be obtained from 

N-2 
n=0 



(0) 
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The precise upper limit (be it N — 2 or e.g. N/2) has 
little importance, because <f>i^ n turn out to decrease 
rapidly for n > n max , where in general n max <C N/2 
(cf. below). 

(iii) According to ref. jT3], the \at \ decrease only up to some 

4nax, meaning that X^=cT l a ^l 2 snows a plateau as a 
function of i? max , but eventually it diverges. Ref. [13] 
suggests choosing £ max from this plateau. We find, 
however, that in practice it is more useful to monitor 



the sum 



e=o 



; the reason is discussed in connection 
with eq. © below, 
(iv) Choosing i max according to some criterion, J + can be 
approximated through 



J+(V) ~ e - e " Y^aeU^e- 1 ) , 

1=0 

where, as usual, the Lg can be constructed from 

L (x) = 1 , Li(x) = 1 - x , 
Li(x) = 2Lt-i(x) - Ln-i(x) 

(1 + x)L t -i(x) - Lt- 2 {x) 



(7) 



I 



In general the asymptotic value, 



J + (oo, ■) = y^ g e 
e=o 



, I > 2 . (8) 



(9) 



does not vanish. In contrast, physically relevant 
current-current correlators should vanish at infinite 
time separation [T3]. It turns out that this supplemen- 
tary information can be used to choose values of l ma x, 
namely those for which J + (oo,-) vanishes approxi- 
mately, offering "windows of opportunity" , in which 
the algorithm appears to perform reasonably well even 
with non-ideal data, 
(v) The spectral function can be obtained as 



P(&,-) 



disin(wf) J + (i, -)~J + (oo r ) , (10) 



where we have subtracted by hand any possible "rem- 
nant" J + (oo,-). (Otherwise p(u>,-) would diverge as 
~ 1/Q at small frequencies.) 

As the proofs in ref. |13j show, in the limit N — > oo and 
vanishing errors these steps do yield the correct spectral 
function for current-current correlators of the considered 
type. 



thereby renders problem (ii) worse. A finite volume can in 
principle also change the behaviour of Euclidean correla- 
tors in a physically interesting way (cf. ref. [16] and ref- 
erences therein), however we wish to ignore these effects 
for now. (Formally a smooth infinite-volume type shape 
could be obtained e.g. by considering a suitable Gaussian 
smoothing of p(lj, ■)/&.) 

(i) In order for the recipe to apply, we must be able to 
compute the Fourier coefficients G(ui n , •); this implies 
that Q(t, ■) must be integrable around r = mod /3. 
In fact, as mentioned above, in ref. |13] G(t,-) was 
even assumed to be continuous (and therefore finite) 
at r = mod j3. In contrast, the correlation functions 
of composite operators that are relevant for the deter- 
mination of transport coefficients or particle produc- 
tion rates in QCD diverge at small r in the continuum 
limit. 

In principle, a possible way around the problem might 
be to consider the temperature derivative of a cor- 
relator, rather than a correlator as such. This could 
work if the derivative is taken in fixed physical units. 
Sometimes, such differences are rather taken in scaled 
units, i.e. by subtracting values of Q(r, •) = (3 p Q(ff3, •), 
< f < 1, at different (3's, but such differences are 
continuous only if there are no scaling violations at 
small t. 

In terms of a spectral function, the finiteness of G(0 + , ■) 
necessitates p(uj, ■)< C/{to ln 2 w) at u T, cf. eq. (fl"2")l . 
The asymptotics of various spectral functions in this 
regime have been analyzed in ref. [17] . and for vector 
current correlators the decay of the thermal part is 
indeed fast enough to satisfy the bound. 

(ii) As can be seen from eq. ©, the factor (—1)™ implies 
a cancellation between Fourier modes in the construc- 
tion of the a/s; when the a/s multiply the Laguerre 
polynomials in eq. (J7J, further cancellations take place, 
particularly at large t (cf. eq. ©). This implies a 
substantial significance loss. Furthermore, as can be 
seen from eq. ([5]), for a fixed t the coefficients 4>e,n 
grow very fast at small n, reaching a maximal value 
at 

^max ~ -\/2(^ + 1). The maximal value, <^.n max ; 
grows with £, exceeding 10 4 at I = 21 (for n max = 6), 
10 6 at I = 40 (for n max = 9) and 10 8 at £ = 64 
(for n max = 11). So, the Fourier coefficients around 
n ~ "-max would need to be determined with the corre- 
sponding relative accuracy in order to obtain a mean- 
ingful signal even after the cancellations. This is obvi- 
ously a formidable challenge, particularly considering 
that problem (i) already requires a subtraction and 
inflicts an associated significance loss. 



5 Problems 

We now turn to two problems that limit the usefulness 
of the algorithm specified above. For simplicity the dis- 
cussion will be carried out in the combined continuum 
and infinite- volume limit; a finite cutoff alleviates problem 
(i) but adds more structure to the spectral function and 



6 Test 

Despite the problems of the previous section, we now wish 
to demonstrate that in principle the method may still work 
on the qualitative level. In order to achieve this, we assume 
that a suitable ultraviolet subtraction has been carried 
out, and that our data are quite precise. 
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Fig. 2. Top: 2~_^™o x Q>1 as a function of f ma x, for TV = 24 (left) and N = 48 (right); here = ae/T 3 and <r indicates the local 
relative standard deviation of Q(t, •) (results are shown for one random configuration). Bottom: the corresponding X^™(T \® e \ 2 - 



As a model, we take inspiration from a spectral func- 
tion related to an "electric field" correlator yielding the 
momentum-diffusion coefficient of a heavy quark [T8,19, 
20 . We assume a substantial positive intercept of p(ui, -)/uj 
at zero frequency (cf. ref. [21]) and decreasing behaviour 
at large frequency, just fast enough to yield a continuous 



Ptest 



CCo 



(2 + u 2 f 



1 



Co 2 



\n 2 (2 + u 2 ) 



2nT 



(11) 



This spectral function is not positive-definite in order 
to reflect the fact that a suitable ultraviolet subtrac- 
tion has been carried out, and because a negative p ~ 
— T 4 / (w In 2 w) is precisely the qualitative asymptotic be- 
haviour found in perturbation theory |22| (the logarithm 
squared assumes that the gauge coupling is let to run with 
w)0The coefficient C appears linearly in all steps so that, 
without loss of generality, we set C = 4T 3 in the following, 
thereby normalizing p tcs t / (ujT 3 ) to unity at ui — > 0. 

The Euclidean correlator is subsequently integrated 
numerically from 



dcu 



cosh 



"Ptcst" 



(«-')' 



sinh 4p 



(12) 



for r = k/3/N < (3/2. At each r we add a random error 
from a Gaussian distribution of relative variance a 2 , and 
then mirror Q(t, •) to the whole interval through the sym- 
metry in t — y [3 — t. On this "data" the steps of section 4 
are applied. Small variants of eq. (jllj) bring along little 
change, but the situation deteriorates rapidly if the struc- 
ture is more peaked either in the ultraviolet (u> 3> 27rT) 
or in the infrared (w <C 2irT). 



3 Note that the Lorentzian form p ~ u/(rj 2 
decrease fast enough at large frequencies. 



The behaviour of the coefficients at, in particular the 
"diagnostic" sum of eq. ©, is illustrated in fig. [5J the 
quality of recovering the spectral function using windows 
of opportunity deduced from fig.[2]is shown in fig. [3] Fig. [4] 
demonstrates the effect of statistical noise on the final 
results. The lessons we draw are the following: 

(i) The recovery is reasonable only when X^=o* a e ~ 0- 
For good accuracy, the sum X^=cT a £ snows a near- 
zero minimum; then ^ max should be chosen close to 
the minimum (the function p(u>, is also extremal 
there). For poor accuracy it could happen that no clear 
cf. a = 10" 3 in fig. deleft); then 



minimum is seen, 
f ra should be chosen close to a point where ^2 



would not 



In 

1=0 



crosses zero. 

(ii) Other values of € max lead in general to nonsensical re- 
sults. 

(iii) Within a given "robust" near-zero minimum, the de- 
pendence on N and a is quite mild. 

(iv) The recovery can be qualitatively improved only by 
increasing the accuracy so much that a second window 
opens up (cf. a — 10 -8 in fig.H]). This is unlikely to be 
reached in practice and, in any case, the improvement 
is not that overwhelming (cf. fig. [3]). 

(v) The statements above apply to any single random con- 
figuration. With a sample of them, £ max could be sep- 
arately fixed for each configuration. Within our toy 
model, it requires an accuracy a < 10 to find a use- 
ful minimum for almost every configuration (cf. fig-H]). 
For a — 10~ 3 , typical configurations show no near-zero 
minimum (cf. fig. [5]), but rare ones do and if it is pos- 
sible to restrict the statistics to those and to carry out 
the averaging on the final p(uj,-)/w, then a rough esti- 
mate can still be obtained. 

Although it goes beyond the scope of the present note 
to carry out a detailed investigation of issues related to 
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Fig. 3. ptcst/i* = Ptest/(i , r 3 ) as a function of uj = uj/(2ttT), for various relative accuracies a, as well as ^ max chosen from 
"windows of opportunity" in which X^=(T ® e approximately vanishes, cf. fig. [2] (results are shown for one random configuration). 
The thick solid line is the correct (input) result. The case a 
in X^™q x ® e but ^ i s n °t near zero (cf. fig. 0. 



10 with N = 48 shows a "failed" example: there is a minimum 



statistical analysis, we note that, in general, the output 
function depends non-linearly on input data, because the 
value of ^ max varies and affects significantly the result. Er- 
ror estimation should therefore be carried out with e.g. 
jackknifc or bootstrap methods, perhaps with blocked 
configurations (the effect of blocking has been shown to 
be beneficial in connection with the Maximum Entropy 
Method, see e.g. ref. g3]). 

7 Conclusions 

The algorithm of ref. [T3] possesses a number of attractive 
features: it can be fully specified in a small number of ex- 
plicit steps; it requires no priors; it does not necessitate a 
positive-definite spectral function; and it projects out the 
Matsubara zero-mode contribution whose handling has 
been considered a problem in certain contexts. 

Unfortunately, from a practical point of view, the 
algorithm of ref. [13 cannot be guaranteed to yield a 
quantitatively accurate analytic continuation of thermal 
imaginary-time data. In some sense, the situation is akin 
to the sign problem hampering simulations of QCD with a 
finite baryon number density: there are significant cancel- 
lations taking place, particularly if a spectral function at 
a small frequency u> <C 2irT needs to be determined. Also, 
short-distance divergences need to be subtracted from the 
Euclidean correlator Q{t, •), which constitutes a signifi- 
cance loss of its own. 

Nevertheless, we have demonstrated that in a lucky 
case with a structureless spectral function and precise data 
(with relative errors < 0.1% after the ultraviolet subtrac- 
tion), already N > 20 data points may yield a qualitative 
reproduction of a transport coefficient (zero-frequency in- 
tercept of p(uj, In general, it is difficult to estimate 



systematic errors, but if a clear near-zero minimum in 
Y^i=q a i i s found as a function of € max , then it appears 
that a < 50% uncertainty can be expected. This could al- 
ready be useful, given that current model-independent de- 
terminations of transport coefficients might contain errors 
of more than 100% [2T] . 
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